//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry                         : env klmperry

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"

// abc
cd $abcjpeanalysis
use append-abccare_iv.dta, clear
keep if program == "abc"
egen iq0  = rowmean(sb4y sb5y)
egen iq5 = rowmean(sb5y iq12y iq8y iq15y)
egen iq2 = rowmean(iq6y6m iq7y iq6y)
rename iq21y iq15
keep id treat iq0 iq2 iq5 iq15 male
gen abc = 1
tempfile abc
save "`abc'", replace

// perry
cd  $perrydatas
use perry_base.dta, clear
gen     treat = 1 if C2 == 2
replace treat = 0 if C2 == 1
gen     male = 1 if C4 == 1
replace male = 0 if C4 == 2 
gen perry = 1
keep    id treat sb_iq_5 sb_iq_7 sb_iq_10 perry male
rename sb_iq_5  iq0 
rename sb_iq_7  iq2 
rename sb_iq_10 iq5
gen    iq15 =0
append using "`abc'"

global sex0 & male == 0 
global sex1 & male == 1
global sex2
foreach gen in 0 1 2 {
	foreach prog in perry abc {
		foreach num of numlist 0 1 {
			summ   iq0 if treat == `num' & `prog' == 1 ${sex`gen'}
			matrix sd`prog'_`num'_g`gen' = r(sd) 
		}
		reg iq0 treat if `prog' == 1 ${sex`gen'}
		matrix md`prog'_g`gen' = e(b)
		matrix md`prog'_g`gen' = md`prog'_g`gen'[1,1]
	}
}

// this is an example with the data you constructed for me
cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

global inputs_std cum_avg_daycare_36m_sum parenting_ages13_std
keep   ihdp site tg bwg $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs

// power calculations for singletons and twins by gender, based on perry parameters
gen ones = 1
rename sex male 
levelsof site, local(sites)
foreach sib in 0 1 {
	foreach gen in 0 1 2 {
		local nsites  = 0
		local m1sites = 0 
		local m2sites = 0 
		foreach  site in `sites' {
			local nsites = `nsites' + 1
			summ  ones if twin == `sib' & tg == 0 & site == "`site'" ${sex`gen'}
			local m1sites = `m1sites' + r(sum) 
			summ  ones if twin == `sib' & tg == 1 & site == "`site'" ${sex`gen'}
			local m2sites = `m2sites' + r(sum) 
		}
		foreach m in 1 2 {
			local m`m'sites = round(`m`m'sites'/`nsites',1)
		}
		// use perry parameters
		foreach prog in abc perry {
			local sd`prog'_0 = sd`prog'_`sib'_g`gen'[1,1]
			power twomeans 0, sd(`sd`prog'_0') m1(`m1sites') m2(`m2sites') k1(`nsites') k2(`nsites') rho(0.03) power(0.55(0.05).95) graph 
			matrix power_s`prog'_s`sib'_g`gen' = r(pss_table)
			matrix power_s`prog'_s`sib'_g`gen' = [power_s`prog'_s`sib'_g`gen'[1...,2],power_s`prog'_s`sib'_g`gen'[1...,7]]
			matrix colnames power_s`prog'_s`sib'_g`gen' = po_s`prog'_s`sib'_g`gen' de_s`prog'_s`sib'_g`gen'
			svmat  power_s`prog'_s`sib'_g`gen', names(col)
		}
	}
}

// graph
keep po_sabc_s0_g0-de_sperry_s1_g2
drop if po_sabc_s0_g0 ==.
// set treatment effects of perry 
foreach gen of numlist 0 1 2 {
	gen mdperry_g`gen' = mdperry_g`gen'[1,1]
	gen mdabc_g`gen'   = mdabc_g`gen'[1,1]
}


// pooled
cd $output
#delimit
twoway         (scatter  de_sperry_s0_g2 po_sperry_s0_g2, msize(medium)  mcolor(gs0)   msymbol(circle)  connect(l) lcolor(gs0)   lwidth(medthick) yline(2 , lpattern(solid) lcolor(none)))
	       (scatter  de_sperry_s1_g2 po_sperry_s0_g2, msize(medium)  mcolor(gs10)  msymbol(diamond) connect(l) lcolor(gs10)  lwidth(medthick) yline(10, lpattern(solid) lcolor(none)))
	       
	       (line       mdperry_g2 po_sperry_s0_g2, lcolor(gs5) lwidth(medthick) lpattern(dash) xline(.8, lpattern(dash)  lcolor(gs12) lwidth(vvvthick)))
	       ,

		legend(label(1 "Singletons")        label(2 "Twins")        
		       label(3 "Perry")
		       
			position(6) size(medsmall) order(1 2 3) rows(1) region(lcolor(gs0)))
		  xlabel(.55[.1].95 .8, labsize(medium)        grid   glcolor(none) glpattern(solid)) 
		  ylabel(2.5[2]10.5,   labsize(medium) angle(h)     glcolor(none) glpattern(solid)) 
		  ytitle("Treatment Effect", size(medium))
		  xtitle("Power")
		  graphregion(color(white)) plotregion(color(white));
#delimit cr
graph export powerall_perry.eps, replace

#delimit
twoway         (scatter  de_sabc_s0_g2 po_sabc_s0_g2, msize(medium)  mcolor(gs0)   msymbol(circle)  connect(l) lcolor(gs0)   lwidth(medthick) yline(2 , lpattern(solid) lcolor(none)))
	       (scatter  de_sabc_s1_g2 po_sabc_s0_g2, msize(medium)  mcolor(gs10)  msymbol(diamond) connect(l) lcolor(gs10)  lwidth(medthick) yline(10, lpattern(solid) lcolor(none)))
	       
	       (line       mdabc_g2 po_sabc_s0_g2, lcolor(gs5) lwidth(medthick) lpattern(dash) xline(.8, lpattern(dash)  lcolor(gs12) lwidth(vvvthick)))
	       ,

		legend(label(1 "Singletons")        label(2 "Twins")        
		       label(3 "ABC")
		       
			position(6) size(medsmall) order(1 2 3) rows(1) region(lcolor(gs0)))
		  xlabel(.55[.1].95 .8, labsize(medium)        grid   glcolor(none) glpattern(solid)) 
		  ylabel(2.5[2]10.5,   labsize(medium) angle(h)     glcolor(none) glpattern(solid)) 
		  ytitle("Treatment Effect", size(medium))
		  xtitle("Power")
		  graphregion(color(white)) plotregion(color(white));
#delimit cr
graph export powerall_abc.eps, replace

// by gender
foreach sib of numlist 1 0 {
	#delimit
	twoway         (scatter  de_sperry_s`sib'_g1 po_sperry_s`sib'_g1, msize(medium)  mcolor(gs0)  msymbol(circle)  connect(l) lcolor(gs0)  lwidth(medthick) yline(2 , lpattern(solid) lcolor(none)))
		       (scatter  de_sperry_s`sib'_g0 po_sperry_s`sib'_g0, msize(medium)  mcolor(gs5)  msymbol(diamond) connect(l) lcolor(gs5)  lwidth(medthick) yline(10, lpattern(solid) lcolor(none)))
		       
		       (line     mdperry_g1 po_sperry_s0_g1, lcolor(gs0)  lwidth(medthick) lpattern(dash) xline(.8, lpattern(dash)  lcolor(gs12) lwidth(vvvthick)))
		       (line     mdperry_g0 po_sperry_s0_g0, lcolor(gs5)  lwidth(medthick) lpattern(dash)),
		      
			legend(label(1 "Male")        label(2 "Female")        
			       label(3 "Perry: Male") label(4 "Perry: Female")
			       
				position(6) size(medsmall) order(1 2 3 4) rows(2) region(lcolor(gs0)))
			  xlabel(.55[.1].95 .8, labsize(medium)        grid   glcolor(none) glpattern(solid)) 
			  ylabel(2[4]18,   labsize(medium) angle(h)     glcolor(none) glpattern(solid)) 
			  ytitle("Treatment Effect", size(medium))
			  xtitle("Power")
			  graphregion(color(white)) plotregion(color(white));
	#delimit cr
	graph export powerbygender_s`sib'.eps, replace
}
